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SUMMARY 


Supercritical wing technology is expected to have a signifi- 
cant influence on the next generation of commercial aircraft. 
Computational fluid dynamics is playing a central role in the 
development of new supercritical wing sections. One of the prin- 
cipal tools is a fast and reliable code that simulates two- 
dimensional wind tunnel data for transonic flow at high Reynolds 
numbers. This is used widely by industry to assess drag creep and 
drag rise. Codes for the design of shockless airfoils by the 
hodograph method have not been so well received because they 
usually require a lot of trial and error. However, a more advanc- 
ed mathematical approach makes it possible to assign the pressure 
as a function of the arc length and then obtain a shockless air- 
foil that nearly achieves the given distribution of pressure. 

This tool should enable engineers to design families of transonic 
airfoils more easily both for airplane wings and for compressor 
blades in cascade. 


INTRODUCTION 


There are plans to use the supercritical wing on the next 
generation of commercial aircraft so as to economize on fuel 
consumption by reducing drag. Computer codes have served well in 
meeting the consequent demand for new wing sections. One of the 
most widely adopted codes was developed at the Courant Institute 
to simulate two-dimensional transonic flow over an airfoil at high 
Reynolds numbers (ref. 1) . This work is an example of the possibil- 
ity of replacing wind tunnel tests by computational fluid 
dynamics. Another approach to the supercritical wing is through 
shockless airfoils. Here a novel boundary value problem in the 
hodograph plane will be discussed that enables one to design a 
shockless airfoil so that its pressure distribution very nearly 
takes on data that have been prescribed. An advanced design code 
of this kind has been written recently by David Korn and is turning 
out to be so successful that it may ultimately gain the same 
acceptance as the better established analysis code. 
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Physically realistic transonic flow computations can be based 
on a potential equation that presupposes conservation of entropy 
across shock waves, but permits a jump in the normal component of 
momentum. However, to treat either the problem of design or of 
analysis for transonic airfoils in a satisfactory way from the 
engineering point of view, it is necessary to take into account 
the effect of the turbulent boundary layer. The simplest proce- 
dure is to calculate the displacement thickness of the boundary 
layer from the inviscid pressure distribution by a momentum inte- 
gral method of Nash and Macdonald (ref. 2) . For analysis one adds 
the displacement thickness to the profile at each cycle of an 
iterative scheme determining the flow. In the case of design a 
corresponding quantity is subtracted from the airfoil coordinates, 
which therefore have to be provided with a slightly open trailing 
edge to begin with. 

It is important to eliminate separation entirely in the prob- 
lem of design if there is to be no loss of lift in practice. This 
can be accomplished by imposing a pressure distribution at the 
rear of the upper surface that just avoids separation according to 
a criterion of Stratford (ref. 3) . The boundary layer correction 
has been found to give satisfactory results even when its imple- 
mentation only involves a primitive model of the wake in which 
pressure forces balance across a parallel pair of trailing stream- 
lines. Extensive wind tunnel tests from laboratories all over 
the world confirm that the analysis code agrees well with experi- 
mental data when the boundary layer correction is made. Prelimin- 
ary test data on a cascade airfoil that was heavily aft-loaded 
also inspire confidence in the concept of using a Stratford pres- 
sure distribution to avoid loss of lift in design by the hodograph 
method. 

The transonic flow codes developed at the Courant Institute 
have been distributed to industry by the Langley Research Center. 

In the future they will also become available through the Argonne 
Code Center of the Argonne National Laboratory. 

THE METHOD OF COMPLEX CHARACTERISTICS 


The partial differential equations for the velocity potential 
<j> and stream function ^'of two-dimensional irrotational flow of a 
compressible fluid can be written in terms of characteristic 
coordinates £ and ri in the canonical form 

<p£ = i i/i-M 2 i/^/p , <f> n = -i/l-M 2 ^ /p , 

where the local Mach number M and the density p are functions of 
the speed q defined by Bernoulli's law. A fast and flexible 
numerical scheme for the construction of smooth transonic flows in 
the hodograph plane has been developed by continuing these equa- 
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tions analytically into the domain of complex values of the two 
independent variables £ and n (ref. 4) . The coordinates % and n 
can be specified in terms of the speed q and the flow angle 9 by 
the formulas 

log f(£) = | /l-M 2 5a - i9 , log f(n) = J /l-M 2 ^SL + i6 , 

where f is any complex analytic function. Prescription of a sec- 
ond arbitrary analytic function g serves to determine <J> and i/> as 
solutions of the characteristic initial value problem 


iM£rn 0 ) = g(£) , «MS 0 *n) = g(n) , 

where = ^ is a fixed point in the complex plane. With these 

conventions it turns out that tJj(£,ti) - iMh*!)/ whence for sub- 
sonic flow the real hodograph plane corresponds to points in the 
complex domain where £ = h . 

Consider the nonlinear boundary value problem of designing 
an airfoil on which the speed q has been assigned as a function 
of the arc length s. To construct a solution it is helpful to 
view f as a function mapping the region of flow onto the unit 
circle |?| <1. There both log f and g have natural expansions 
as power series in £ after appropriate singularities accounting 
for the flow at infinity have been subtracted off. The coeffici- 
ents of truncations of these series can be determined by inter- 
polating to meet boundary conditions on q and iJj at equally 
spaced points of the circumference |£| = 1. Such a numerical 
solution is easily calculated because the matrix of the system of 
linear equations for the coefficients is well conditioned. This 
analytical procedure has the advantage that its formulation can 
be extended to the case of transonic flow so as to yield a shock- 
less airfoil nearly fitting the prescribed data even when an 
exact solution, of the physical problem does not exist. 

To calculate transonic flows by the method that has been 
proposed, it is necessary to circumvent the sonic locus M = 1, 
which becomes a singularity of the partial differential equations 
for <J> and i/j in canonical form. In the plane £ = n this locus 
separates the region of subsonic flow from a domain where 
is no longer real. In the latter domain it is necessary to 
extend in some empirical fashion the relationship between <J> and 
s that is imposed by assigning q as a function of s. A formula- 
tion of the boundary conditions that applies to both the subsonic 
and the supersonic flow regimes is given by the formulas 

Re {log f(£)> = h , Re {<{/(£,£) > + k Im > - 0 

on | £ | = 1 , where k is a real constant and h is a function of 
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Re {<f>(£jl)} obtained from the known relationships among s, q and 
<j>. The nonlinearity of the problem makes it necessary to iterate 
on this relationship in finding a numerical solution. 

Empirical data on the condition number of the matrix for the 
linear equations determining the power series coefficients of the 
analytic function g indicate that the boundary value problem for 
rp that has been formulated is well posed even in the transonic 
case. In contrast with the Tricomi problem, boundary values are 
assigned around the whole circumference of the unit circle. The 
success of the procedure can be attributed to the fact that data 
are assigned in a suitable complex extension of the real plane. 

In general limiting lines may appear in the physical plane, 
but it has been found that these can be suppressed by appropriate 
selection of the rules defining the function h and the real para- 
meter k that occur in the specification of the boundary condi- 
tions. Thus a tool becomes available for the construction of 
supercritical wing sections from their pressure distributions. 
Figure 1 shows an example of a shockless airfoil that was obtained 
this way, together with its Mach lines. Observe that the input 
pressure coefficient Cp differs somewhat from the values 
calculated as output of the flow in the supersonic zone, which is 
rather large. The data that were assigned are based on a modifi- 
cation of the experimental pressure distribution on Whitcomb's 
original supercritical wing (ref. 5) shown in Figure 2. 

The design code has been written to include the case of tran- 
sonic airfoils in cascade. This model seems to offer consider- 
able promise for improvement in the efficiency of certain stages 
of high speed compressors. However, to handle cascades of high 
solidity with adequate resolution it is desirable to replace a 
conformal mapping onto the unit circle [ £ | < 1 by the mapping 
onto an ellipse, where the Tchebycheff polynomials become prefer- 
able to powers of £ for expansion of the analytic functions 
log f and g. Likewise, to achieve adequate resolution at the 
trailing edge in cases of heavy aft-loading it is helpful to 
insert a special term at the tail in the representation of the map 
function f. 

The new code represents a major advance over what was 
achieved in earlier versions, whose use required excessive trial 
and error (ref. 4) . A typical run takes about six minutes on the 
CDC 6600 computer. Closure of the airfoil is readily attained by 
adjusting the pressure at the trailing edge and the relative 
lengths of arc over the upper and lower surfaces between the 
stagnation point and the trailing edge. A general principle to 
be observed when using shockless airfoils to design supercritical 
wing sections is that drag creep can be reduced by diminishing the 
size of the supersonic zone of flow, especially toward the rear of 
the profile. In practice the best way to assess the performance 
of a new design is to run it through the analysis code, which will 
be discussed next. 
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ESTIMATION OF THE DRAG 


Analysis of the transonic flow past an airfoil can be based 
on a partial differential equation for the velocity potential <j>. 
Weak solutions modelling shock waves are calculated by adding 
artificial viscosity. This can be accomplished with a full con- 
servation form (FCF) of the equation, but a simpler version not 
in conservation form (NCF) is sometimes more useful (ref. 1) . To 
handle the boundary conditions it is convenient to map the region 
of flow conformally onto the exterior of the unit circle. If r 
and 0 stand for polar coordinates there, the quasilinear equation 
for <f> can be written as 


a <f >00 + 2 b <j> 0r + c <f> rr + d = 0 

when artificial viscosity is omitted. The simplest way of 
introducing artificial viscosity numerically, suggested first by 
Murman and Cole in a fundamental paper (ref. 6 ) , is to use finite 
difference approximations that are retarded in the direction of 
the flow, which for practical purposes can be taken as the direc- 
tion of 6 . This does not perturb the Neumann boundary condition 
on cf) . 

The finite difference equations for transonic flow can be 
solved iteratively by a variety of relaxation schemes, all ojl 
which take the form of marching processes with respect to an arti- 
ficial time parameter. Antony Jameson has found that the rate of 
convergence can be accelerated by substituting a fast solver over 
the subsonic flow region between every few cycles of relaxation 
(ref. 7) . Such a procedure has been programmed by Frances Bauer 
using fast Fourier transform with respect to the periodic variable 
0. This reduces the calculation time by a factor of three even 
when a boundary layer correction is included in the computation. 

A standard run of her airfoil code now takes less than three min- 
utes on the CDC 6600 computer. 

Detailed comparisons with experimental data show that the 
NCF transonic equation gives significantly better simulation of 
shock wave-boundary layer interaction than does the FCF equation, 
especially in cases with a shock at the rear of the profile where 
the turbulent boundary layer is relatively thick. It would appear 
that the NCF method leads to less radical gradients in the pres- 
sure behind the shock, which is consistent with the observations. 
The NCF and experimental speeds both tend to jump down barely 
below the speed of sound behind a shock. Figure 2 shows the kind 
of agreement between theoretical and test data that is usually 
seen. Wall effect is accounted for by running the computer code 
at the same lift coefficient C T that occurs in the experiment. 

JLi 

Because of erroneous positive terms in the artificial viscos- 
ity, the shock jumps defined by the NCF method create mass instead 
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of conserving it. However, the amount of mass produced is only of 
the order of magnitude of the square of the shock strength for 
nearly sonic flow. The resulting errors are therefore negligible 
except for their effect on the calculation of the wave drag, which 
has the order of magnitude of the cube of the shock strength. A 
correct estimate of the drag can be obtained from NCF computations 
by working with the path- independent momentum integral 

D = J [p dy + (<f> x - cjdifi] , 


where p and c* stand for the pressure and the critical speed, 
respectively. The integrand has been arranged so that across a 
normal shock wave parallel to the y-axis it jumps by an amount of 
the third order in the shock strength. Therefore integration 
around the shocks gives a reasonable measure of the wave drag 
even when mass is not conserved. 


The path of integration can be deformed onto the profile to 
define a standard integral of the pressure there, but a correction 
term evaluated over a large circle should be added because of a 
sink at infinity accounting for the mass generated by the NCF 
method. Let £, p and denote the chord length of the airfoil, 

the density at infinity and the speed at infinity, respectively. 
The corrected formula for the wave drag coefficient C DW becomes 


' DW * P q 2 


p dy - 2 


c *-q 0 

* 


I 


d^ , 


where the first integral is extended over the profile and the 
second integral is extended ovef a large circle separating the 
profile from infinity. In Figure 3 a comparison is presented 
between experimental, corrected NCF, uncorrected NCF and FCF 
values of the total drag coefficient C_ for a shockless airfoil 
tested at Reynolds number R = 20><10 6 by Jerzy Kacprzinski at the 
National Aeronautical Establishment in Ottawa. The corrected drag 
formula is seen to give a fairly reliable assessment of the per- 
formance of the airfoil. 


There are examples where the results of the NCF code agree 
well with experimental data right up to the onset of buffet. 

Shock locations are predicted with remarkable accuracy over a wide 
range of conditions, although some improvement would be desirable 
at lower Reynolds numbers where transition becomes important. Thus 
the analysis code has been adequately validated for simulation of 
experimental data in two-dimensional flow. In particular, it 
models the trailing edge in a satisfactory way even for heavily 
aft-loaded airfoils. It is therefore of some interest that the 
code predicts no loss of lift for airfoils designed by the hodo- 
graph method when a Stratford distribution is used to eliminate 
separation completely over the whole profile. It would neverthe- 
less be desirable to confirm this result experimentally by further 
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testing of shockless airfoils such as the one shown in Figure 1. 

There is need for more research on computational methods for 
transonic flow. The progress in supercritical wing technology 
should be extended to cascades of airfoils and flows in turbo- 
machinery. For the immediate future, the most challenging problem 
is analysis of the flow past wing-body combinations modelling an 
airplane in three dimensions. As a first step it would seem that 
the NCF equation for a velocity potential furnishes the most 
feasible mathematical formulation. Perhaps the Bateman variation- 
al principle asserting that the volume integral of the pressure is 
a stationary functional of the velocity potential, applied in the 
context of the finite element method, offers the best prospect of 
deriving convenient difference equations, provided artificial 
viscosity can be added successfully. 
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